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Abstract 

Actively propelled particles undergoing dissipative collisions are known to develop a state of 
spatially distributed coherently moving clusters. For densities larger than a characteristic value 
clusters grow in time and form a stationary well-ordered state of coherent macroscopic motion. In 
this work we address two questions: (i) What is the role of the particles' aspect ratio in the context 
of cluster formation, and does the particle shape affect the system's behavior on hydrodynamic 
scales? (ii) To what extent does particle conservation influence pattern formation? To answer these 
questions we suggest a simple kinetic model permitting to depict some of the interaction properties 
between freely moving particles and particles integrated in clusters. To this end, we introduce two 
particle species: single and cluster particles. Specifically, we account for coalescence of clusters 
from single particles, assembly of single particles on existing clusters, collisions between clusters, 
and cluster disassembly. Coarse-graining our kinetic model, (i) we demonstrate that particle shape 
(i.e. aspect ratio) shifts the scale of the transition density, but does not impact the instabilities 
at the ordering threshold, (ii) We show that the validity of particle conservation determines the 
existence of a longitudinal instability, which tends to amplify density heterogeneities locally, and 
in turn triggers a wave pattern with wave vectors parallel to the axis of macroscopic order. If the 
system is in contact with a particle reservoir this instability vanishes due to a compensation of 
density heterogeneities. 
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I. INTRODUCTION 



The emergence of collective motion is a ubiquitous phenomenon in nature, encountered 
in a great variety of actively propelled systems [IHS]. Coherently moving groups have been 
observed over a broad range of length scales, spanning from micrometer-sized systems 
[TO] over millimeter large granules [TlHT3] to large groups of animals [14j. The fact that the 
capability of synchronizing movements between agents is shared even among fundamentally 
different systems has called for abstract modeling approaches, aiming at identifying the 
essential properties of these systems both, in terms of analytical descriptions [T3H2H], and 
by means of agent based simulation techniques [2SH3S]. 

Theoretically, the emergence of collective motion has mostly been studied in the context of 
particle conserving systems. There are, however, a number of experimental systems, in which 
the assumption of particle conservation is questionable. In typical gliding assays [lH6l [U [10] , 
for instance, collective motion of filaments is observed on a two-dimensional "motor carpet" 
which itself is in contact with a three-dimensional bulk reservoir of filaments. However, the 
impact of particle conservation on the formation of patterns of collective motion remains 
largely elusive. 

Here, we address the significance of constraints for particle number by highlighting the 
differences in the collective properties between particle conserving systems and those in 
contact with a particle reservoir. Our focus will be on the comparison of two archetypi- 
cal scenarios, which we will refer to as the canonical (particle conserving), and the grand 
canonical (violating particle conservation) scenario, respectively. 

To this end, we will resort to a kinetic approach, which has been set up previously by 
Aranson et al. [16] to describe pattern formation in a system of interacting microtubules, and 
which has been extended to the case of self-propelled spheres by Bertin et al. [T71[22]. In the 
following we will extend this description in accordance with a physical picture of collective 
motion that has been developed over the last decade based on observations in agent-based 
simulations of locally interacting, particle conserving systems |3ll |32l ISTJ [39] . Among the 
most pertinent phenomena that have been reported in the context of these studies is the 
formation of intricate local structures pervading these systems in the vicinity of the ordering 
transition: Densely packed cohorts of coherently moving particles — subsequently referred to 
as clusters — incessantly "nucleate" and "evaporate" on local scales, even below threshold. 
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FIG. 1: Illustration of the canonical and grand canonical modeling framework, highlighting the 
quintessential differences in the context of pattern formation. In the homogeneously polarized 
state (left), the cluster particles density (blue arrows) constitutes the system's macroscopic net 
momentum go, while some fraction of the system's particles, the single particles (orange dots) 
exhibit zero net momentum. Spatial perturbations of both density fields lead to two fundamen- 
tally different outcomes: (i) In case of a closed system obeying total particle conservation (single 
particles + cluster particles), termed as the canonical model, the homogeneously polarized state is 
longitudinally unstable, with a wave vector q parallel to the polarized state go, potentially enforc- 
ing a wave-like pattern, (ii) In contrast, open systems turn out to be stable against this kind of 
density fluctuations. 

rendering the system isotropic and homogeneous only in the limit of macroscopic length 
scales. Individual particles exhibit superdiffusive behavior in this regime, performing quasi- 
ballistic "flights" as long as they are part of a cluster, and conventional particle diffusion if 
they are not. Above threshold, collective motion manifests itself on macroscopic scales in 
the form of coherently moving and dense bands, which are submersed in an isotropic low- 
density "particle sea". Spatially homogeneous flowing states, in contrast, are observed only 
well beyond the ordering threshold [31]. Moreover, particle geometry was demonstrated 
to play an essential role in the context of clustering dynamics, with higher aspect ratios 
facilitating the formation of clusters of coherently moving particles [37] . 

In the light of the above, we suggest a simplified modeling framework to incorporate the 
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intricate role of clusters on the ordering behavior, which will be presented in greater technical 
detail in the following section: Particles interact via binary collisions with a scattering cross 
section which is explicitly derived as a function of particle shape. Depending on whether a 
given particle is part of a cluster or not, it will be associated with one of two distinct particle 
classes, which we will refer to as the class of cluster particles and the class of single particles, 
respectively. Single particles are "converted" to cluster particles by "condensation" every 
time a single particle collides with a cluster. Conversely, cluster particles are "converted" 
back to single particles by an "evaporation" process which we assume to occur at some 
constant (possibly particle shape dependent rate. Moreover, in the absence of inter- 
actions, cluster particles will be assumed to move ballistically, whereas single particles will 
be assumed to perform random walks. Taken together, the conversion dynamics and the 
class-specificity of particle motion provide a simple way to implement the typical superdif- 
fusive behavior of individual particles, that was alluded to above. To assess the importance 
of particle conservation in the context of pattern formation, we will analyze two variants of 
this model: Firstly, we study closed systems in which the total number of particles is con- 
served {canonical scenario) and where, consequently, the denser cluster phase grows at the 
expense of the single phase. Secondly, we examine open systems in contact with a particle 
reservoir {grand canonical scenario), where the particle current out of the single phase is 
compensated so as to retain the density of the isotropic sea of single particles at a constant 
level; cf. figure [TJ 

Our work is structured as follows: In |Tl] the modeling framework for the canonical and 
grand canonical model is introduced and the model equations are discussed in detail. The 
corresponding hydrodynamic equations are derived in |III| by means of an appropriate trun- 
cation scheme in Fourier space. Therein, we also give explicit expressions of the kinetic 
coefficients as a function of the particles' aspect ratio and velocity, noise level and density 



for single particles and cluster particles. IV is devoted to the analysis of the homogeneous 
equations. The dynamic's stationary fixed points are determined and the phase boundary 
between the isotropic and homogeneous state is calculated. |V] deals with the implications 
of the inhomogeneous equations in the framework of a linear stability analysis, which are 
concluded in IVII 
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II. COARSE-GRAINED KINETIC MODEL 



We consider rod-like particles of length L and diameter d moving in two dimensions with 
a constant velocity v. A particle's state is determined by its position x and the orientation 
9 of its velocity vector. To describe the time evolution of the system, we adopt a kinetic 
approach [I6HIH1E2]. 

On mesoscopic scales, the system's spatio-temporal evolution is then governed by 
Boltzmann-like equations for the one-particle distribution functions within the classes of 
single particles and cluster particles, respectively. Interactions enter this description by 
means of collision integrals. The kernel of these integrals involves both, a measure for the 
rate of collisions, as well as a "collision rule" implementing a mapping between pre- and 
post-coUisional directions 6 and 6' of each of the two partaking particles. Here, we are led to 
consider a simplified model of binary particle interactions, which builds on the distinction 
between single particles and cluster particles. The details of this model will be described in 
the following section. 

A. Reaction equations 

Let S{6) and C{6) refer to a particle moving in the direction of 6 and being associated with 
the class of single particles or cluster particles, respectively. In the absence of interactions, 
single particles are assumed to perform a persistent random walk, which we model as a 
succession of ballistic straight flights, interspersed by self-diffusion ("tumble") events. These 
tumble events are assumed to occur at a constant rate A and reorient the particle's orientation 
by a random amount t^o- 



with (To denoting the standard deviation. On times scales much larger than A^^, this tum- 
bling behavior can be described as conventional particle diffusion, with the particles' diffusion 
constant being a function of A and ctq [40J . 




(1) 



For simplicity we assume -i^o to be Gaussian-distributed, 




(2) 
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FIG. 2: (a) Illustration of two single particle species (light orange) with a pre-collisional relative 
angle of 612, colliding such that they align collinear to the average angle 9. Both particles become 
a cluster species after the collision, (b) Right: Illustration of a possible scenario, where a single 
particle joins a cluster by perfectly aligning to the cluster particles (blue). Left: A particle leaves 
the cluster by a random change of its direction at a characteristic rate e. 

When two single particles S{6i) and 3(62) collide they are assumed to assemble a cluster, 
i.e. each of the two particles becomes a cluster particle (see figure |2^): 

3(6^) + 8(62) ^cie + ^) + c{e + ^), (3) 

where [12] 

m, 02) = 1(01 +02) (4) 

denotes the average of both pre-collisional angles 0i and 02, and where ^9 is a random variable 
which we, again, assume to be Gaussian-distributed: 

p(^) = ^i=exp(-^V2a2). (5) 

The rate of binary collisions, such as equation ([s]), are determined by a particle-shape de- 
pendent differential scattering cross section, which will be discussed below; see |II B and 

m 

Collisions involving cluster particles are distinct from single particle events. Due to the 
close spatial proximity of particles within each cluster these collisions correspond to many- 
particle interactions. Needless to say, a detailed description of cluster formation and the en- 
suing particle dynamics represents a highly complex matter, requiring explicit consideration 
of such many-particle interactions. For simplicity, we will resort to the following simplified 
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interaction picture: We assume that (binary) collisions between single particles and cluster 
particles lead to a condensation process during which the single particle aligns to the cluster 
particle without changing the direction of the cluster as a whole: 

5(6*1) + C(e2) ^C(02) + C(^2). (6) 

Eq. ^ thus captures the net effect of collisions between single particles and cluster parti- 
cles, during which multiple collisions, involving neighboring particles belonging to the same 
cluster, stabilize the cluster's direction; cf . figure [2] b for an illustration. 

Collisions among cluster particles is an even more intricate process, since they actually 
depend on size and shape of both colliding clusters, and in general involve multi-particle 
interactions. In the framework of a Boltzmann-like description, correlations in the particle 
distribution are neglected and only binary interactions are considered. The frequency of 
interactions are determined by a geometrical construction called the "Boltzmann cylinder" , 
assuming that particle positions are homogeneously distributed on local scales. With regard 
to many-particle interactions during collisions among cluster particles, we thus have to resort 
to some kind of simplified, binary collision picture. Since our kinetic model lacks any direct 
notion of cluster size or shape, we will stick to the assumption that, on average, collisions 
between cluster particles are devoid of any directional bias, leading to the same type of 
collision rule as for single particles: 

c{e,) + cie2) ^cie + ^) + cie + ^). (7) 

Again, constitutes a Gaussian-distributed random variable given in ([s]). Moreover, due to 
external (e.g. thermal background) and internal (e.g. noisy propelling mechanism) noise, 
cluster particles evaporate to become single particles. In analogy to the self- diffusion of 
single particles, we thus introduce a rate [l3j e characterizing the following evaporation 
process: 

C{9) 4 S{e' = 6 + ^o). (8) 

Also in this case, the strength of the angular changes are Gaussian-distributed according 
to ([2]), and for simplicity we use the same standard deviation ao as for the single particles^ 
persistent random walk. As discussed above, cluster particles are strongly caged due to 
their close proximity to neighboring, coUinearly moving particles. Reorientations of cluster 
particles due to noise are therefore strongly counteracted by realigning particle collisions, 
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rendering cluster particles considerably less susceptible to random fluctuations than single 
particles. Hence, we assume 

e < A, (9) 

which is consistent with the observations in agent-based simulations slightly below the or- 
dering transition [31], flnding coherently moving clusters in an unpolarized background of 
randomly moving particles. In this regime individual particles exhibit superdiffusive behav- 
ior, performing quasi-ballistic "flights" as long as they are part of a cluster, and conventional 
particle diffusion if they are not. 

B. Constitutive equations 

Building on the modeling framework deflned above, we now set up a kinetic description for 
the canonical model. We denote by s{9, x, t) and c{9, x, t) the one-particle distribution func- 
tions within the class of single particles and cluster particles, respectively, i.e. s{6, x, t) dO (Px 
gives the number of single particles located in an inflnitesimal region [x, x + rfx] with ori- 
entations in the interval [9,9 + d9] (and likewise for c{9,^,t) d9 d^x). Both one-particle 
distribution functions are subject to convection due to the propelling velocity v of each 
particle. Moreover, local fluctuations in the one-particle distribution functions due to self- 
diffusion and collision events are to be accounted for. We thus arrive at the following set of 
Boltzmann-like equations for the canonical model: 

(9ts(0,x,t) + v- Vs(0,x,t) = s(0,x,t), (10a) 
9tc(^,x,t) + V ■ Vc(e,x,t) = c(^,x,t), (10b) 

where the source terms s(6',x, t) and c(6',x, t) read 

s = A [Vi:^\9)-Vir\9)]+eV[:'\9)-Ci-\9)-A[s,c;9l (11a) 
c = -eVi-\9)+Ci-^\9)+Ci-'\9)+A[c,s;9]-Ci-\9). (lib) 

They give the net number of single particles and cluster particles entering the phase space 
region du = [x, x + c/x] x [9,9 + d9] per unit time and unit area, respectively. The various 
terms correspond to gain [superscript ^"*"''] and loss [superscript'-^-'] of particles by the following 
processes: 
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(i) Self-dijfusion and evaporation. In these cases the source terms are products of the 
corresponding rates and probabihty densities, with 

= f{9) (12) 

denoting the probabihty density for a particular species / to have a certain angle 6*, and 

T^f\0) = {f{0-^o))o (13) 

denoting the transition probability from 9' = 9 — to 9 averaged over all with respect 
to the Gaussian weight (|2|. Note that, here and in the following, the argument of / is 
understood modulo 27r. 

(ii) Collisions within the same class of particles. The collision integrals, representing the 
processes defined in equations ^ and ([T]), are given by standard expressions p^HTSl [22] 

C\-'\9) = l^j dIf{9')f{9")6(^9{9',9") + ^-9^y (14a) 

Cy\9) = j dXf{9')f{9")5{9' -9). (14b) 

Here (...) denotes an average over d G (—00,00) with respect to the Gaussian weig ht (§ 
and the average angle 9 is given in Eq. (|4]). The integral measure 

I dX {...)= r d9' ['^^ d9"r{L,d,\9' -9"\) {...), (15) 

contains the differential scattering cross section 

9' - 9" 



V{L,d,\9' -9"\) = Adv 



sm 



^^^1^ \.H9' - 9") 



(16) 



2 

characterizing the frequency of collisions (i.e. hard-core interactions) between rod-like par- 
ticles. The scattering function F itself carries all information concerning the shape of the 
particles and is a function of the relative orientation of the colliding particles. Reminiscent 
of the Boltzmann scattering cylinder, F can be derived on the basis of purely geometric con- 
siderations assuming that all spatial coordinates within the cylinder are equally probable; 
for details see El 

(iii) Assembly events of a single particle joining a cluster. These events, represented by 
(|6|, occur through binary collisions between single particles and cluster particles and are 
thus represented by an analogous integral expression: 

A[f,g;9] = J dl f{9')g{9")6{9' - 9). (17) 
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III. DERIVATION OF HYDRODYNAMIC EQUATIONS 



In order to reduce our kinetic description to a set of hydrodynamic equations valid on 
large length and time scales, we follow the well-established procedure of Aranson et al. 



and Bertin et al. [HIES], and analyze the angular dependence of equations (10a) and (10b) 
in Fourier space. Due to the 27r-periodicity in 9, the one particle distribution functions can 
be expanded in Fourier series 



,x,t) 
,x,t) 



^ oo 

n=— oo 
^ oo 

— V c„(x,t)e" 



■ind 



irid 



(18a) 
(18b) 



where 



rf^e*"^(^,x,t), 



d9e 



ind , 



X,t). 



(19a) 
(19b) 



s„(x,t) = 

Cn(x,t) = 

J —TV 

Upon identifying o C, e.g. v ^ f e*^ {v = |v|), the zeroth and first Fourier modes are 
directly connected to the hydrodynamic densities ps {single particle density) and pc {cluster 
particle density), and the corresponding current density gs and gc , i.e. 



p,(x,)f:) = So{x,t), 

pc{^,t) = Co(x,t), 

{x,t) = ps{x,t)us{x,t) = vsi{x,t), 

(x,t) = pc(x,t) Uc(x,t) = t;ci(x,t). 



(20a) 
(20b) 
(20c) 
(20d) 



In equations (20c) and (20d), the "=" signs indicate identification of vectors and complex 



numbers. The quantities u^/c denote the velocities of the macroscopic flow flelds established 
by single particles and cluster particles, respectively. Also note that the second Fourier 
components are proportional to the nematic order parameter within the respective class of 



particles (as reflected by the symmetry of e*^^ under 6* — > ^ + vr). Using equations (18a) and 
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(18b), the Boltzmann-like equations (10a) and (10b) transform to 



dtSk + 



dx (sfc+i + Sk-i) - idy (sfc+i - Sfc_ 



(21a) 



Sk- 



dtCk + 



dx (cfc+i + Cfc_i) - idy (cfc+i - Ck-i] 



(21b) 



fCfc -|- ^ ^ -^n,0 ("^nCfc— n CnCfc_n) ~l~ 6^ ^ In,k (■^n-^A;— n ~l~ CnCk- 



where the colhsion integrals In,k are defined as follows: 

1 

In,k = ^l rf0r(L,(i, 101) COS 

27r 



n — 



(22) 



Note, in particular, that Jo,o gives the total scattering cross section. 



A. Truncation scheme 



Equations (21a) and (21b) constitute an infinite set of coupled equations in Fourier space. 



which are fully equivalent to the Boltzmann-like equations (10a) and (10b). To derive a 



closed set of hydrodynamic equations, we need to consider some additional assumptions, 
allowing us to truncate this infinite Fourier space representation. 

Here, our focus will be on virtually isotropic systems in the vicinity of an ordering tran- 
sition breaking rotational symmetry. In this case, deviations of the one-particle distribution 
functions from the constant distribution ~ l/27r are small and contributions from large wave 



numbers in the Fourier series (21a) and (21b) are negligible. We further consider sufficiently 



dilute systems, in which the number of (binary) particle collisions per unit time and area 
[~ (Pc + PsY is much smaller than the corresponding number of single particle diffu- 



sion events [~ \ps\- Together with e ^ A [equation ([9])], stating that disassembly from a 
cluster is strongly hindered by particle caging, allows us to treat single particle diffusion 
as a fast process. The single particle phase thus acts as an isotropic sea of particles where 
particle orientations (but not necessarily particle densities) are equilibrated, and hence the 
net hydrodynamic flow vanishes [u^ = 0]. Finally, from a dimensional analysis of equations 
( pJal ) and ( [l9b| , together with (|20c|) and ( [20dl ), one flnds Ck/pc ~ 0{\\ic\^ /v^). Near the 
onset of order, where |uc|/f ^ 1, we only consider the density (cq) and polarity (ci) of 
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cluster particles, and use the stationary equation for C2 as a closure relation, neglecting all 
contributions from higher order coefficients. 

In summary, we resort to the following truncation scheme, leading to a set of hydrody- 
namic equations, valid near the onset of the ordering transition: 

Sk = 0, V|A;| > 0, (23a) 
Ck = 0, V|A;| > 2. (23b) 



B. Derivation of the hydrodynamic equations 



With the above truncation scheme, (21a) and (21b) reduce to 



dtSo 
dtco 
dtCi 



(24a) 
(24b) 
(24c) 



+ 
+ 



2e "'/^/i,! - /i,o - /o,o ) Co - e + /o,oSo 



2e-'^'/2/2,i - /i,o - ho 



Cl 



C1C2, 



dtC2 



-^[dx + idy] Cl 



(24d) 



+ 
+ 



2e-2-7i 



,0 — -'2,0 



^,0 ) Co - e + Jo nSo 



C2 



'0,0 



^1,0 



ClCi, 



where we used = f^, since f{9) G M (/ G {s,c}), and where 9^e(a) pm.(a)] denote 



the real [imaginary] part of a. Moreover, as can be seen from the definition in (22), the 
collision integrals In,k only depend on the value \n — k/2\, whence only five of the collision 
integrals appearing in the above equations are independent. These integrals as a function 
of the particle's aspect ratio are evaluated and summarized in table |Tj Also note that the 



entire set of equations (24a) - (24d) is independent of the fast single particle diffusion time 



scale A^^ (and, hence, also of the diffusion noise parameter cxo). In our present approach, A 
has only a conceptual meaning in maintaining a well-mixed particle bath within the class of 
single particles. 

For given particle densities, the time scales governing the dynamics of the polar and 
nematic order parameter fields, represented by ci and C2, are given by the linear coefficients 
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in the second line of (24c) and (24d), respectively. As will be detailed in IV B, the onset of 
collective motion is hallmarked by a change in sign of the linear coefficient in (24c), implying 
a diverging time scale for the dynamics of the polarity field. On the other hand the time 
scale for C2 is finite for all densities which implies that the relaxation of the nematic order 



parameter field is fast compared to the polarity field. This allows us to set dtC2 ~ in (24d). 

In the following it will be convenient to write down equations in dimensionless form. To 
this end, we construct the following characteristic scales: Time and space will be measured 
in units of the cluster evaporation time and length scale 

fe = e"^ and 4 = v/e. (25) 

From the cluster evaporation time scale fe and the total scattering cross section Iq^, we can 
construct the characteristic density scale 

Pb = j^. (26) 

The single particle and cluster particle phase constantly exchange particles at rates that are 
determined by cluster evaporation (e) on the one hand [cluster particles — )■ single particles), 
and cluster nucleation due to particle collisions on the other hand {single particles — )■ cluster 
particles), which occur with a rate ~ p/o,o- Therefore the characteristic density scale pf, 
marks the particle density, where both rates balance. In particular, p/p;, = (p^ + Pc)/Pb 
gives the rate of inter-particle collisions relative to cluster evaporation events. Thus, the 
numerical quantity p/pb provides a direct measure expressing the competition between the 
randomizing effects of noise and the order creating effects of particle collisions, hallmarking 
the onset (and maintenance) of collective motion [29]. 
We thus arrive at the following rescaling scheme 



t 






(27a) 


X 


— > 


X ■ ^e. 


(27b) 


Ps/c 


-)■ 


Ps/c ■ Pb, 


(27c) 


g 


-)■ 


. 4 
g- pb—, 

PbTe 


(27d) 






(27e) 



where the characteristic scales for momentum (g) and scattering cross section {In,k) have 
been constructed from those of time, space, and density. In this rescaling the momentum 
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current density is equal to one if the corresponding fluid element with a characteristic density 
Pb, for which cluster evaporation and nucleation balance, is convected with the particle 
velocity v. 



Then, upon eliminating C2 from (24c) as discussed above, and using the relations between 



Fourier modes and hydrodynamic fields (for details seeQ, (20a) - (20d), equations (24a) 



(24d) give rise to the hydrodynamic equations corresponding to the canonical model. In 



rescaled variables they read: 

Pc - {Ps + Pc) Ps 



dtps 
dtPc 
dtg 



V ■ g - Pc + (Ps + Pc) P, 
1^2 



\Vpc + ^V^g 

2 Au2 



(28a) 
(28b) 
(28c) 



+-(g-V)g+^ 

1^2 ^2 



(V 



(g^) 



+ 



p 



^[PcPs]) - -g^a[pc,Ps 



[(V-g)a[pc,p.]-(Vg + Vg*)a[pc,p. 



where 



^[/,^7] = (5/Z/2)V/+(9,Z/2)V^?, 

and where we have introduced the following abbreviations: 





= 1 - (Ps - Pc) + 


(A,0 - 


2e-''"/'h,ij Pc, 


^2 


= 1 - (Ps - Pc) + 


(^2,0 - 


2e-2-'/i,o) Pc, 


P 


— e — ii,o. 








= A,o + hfi - 26^" 


ay 2 J 

-'2,1 




c± 


= -P±2- 







(29) 

(30a) 
(30b) 
(30c) 
(30d) 
(30e) 



Equations (28a) - (28c) capture the evolution of our canonical model system on a hy- 



drodynamic level. More specifically, (28a) and (28b) describe the spatio-temporal evolution 
of the particle densities ps and pc- Since, by the assumptions underlying our model, no 
macroscopic flow of single particles can build up, only the density of cluster particles (pc) 
is subject to convection. This implies that the genuine hydrodynamic momentum field 
g = gc + gs = gc, is carried solely by the subset of cluster particles. Therefore we omit 
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Integral 



Value 



-^0,0 hfl/hfl h,i/hfl hfl/hfl h,i/Io,i 



8dv{2+0 4+g 3 8+7r(g-l) 6-13g 3 7r(l-g)-8 
3tt 5(2+g) 16 2+g 35(2+0 16 2+g 



TABLE I: Summary of relevant collision integrals I^^k ^ a function of the aspect ratio ^ = L/d, 
where L and d denote particle length and diameter, and where v is the particle velocity. The 
quantities In,k/Iofi depend only weakly on the aspect ratio ^. In particular, the signs of In,k/Io,o 
do not change with ^, leaving all our present conclusions made on the basis of the kinetic coefficients 
qualitatively unchanged. 



the subscript c in (28b) and (28c) and denote g = gc- The dynamics of both densities is, 



moreover, driven by source terms, as determined by the reactions discussed in section II A 



The gain and loss parts in these source terms of pc and ps are exactly balanced, such that 
the total density p = pc + Ps is, conserved. As an aside we note that any distinction between 
single particles and cluster particles is a purely conceptual matter. Experimentally, only the 
total density p and the momentum field g are accessible. 



(28c), governing the evolution of the current density g, can be interpreted as a general- 



ization of the Navier-Stokes equation to active systems. The terms on the right hand side of 



(28c) can be given the following interpretation: In the first line, the first two terms account 



for the local dynamics of g. They play a crucial role in establishing and maintaining a state 
of macroscopic flow, as will be detailed below. The Navier-Stokes equation itself, which 
conserves momentum, is devoid of these terms. In formal analogy to the Navier-Stokes 
equation, the density gradient in the first line together with the last term in the second line 
can be interpreted as a pressure gradient. This effective pressure is given by | (^pc + ^ g^^ , 
when neglecting the density-dependence of i>2. The last term in the first line is analogous 



to the shear stress term in the Navier-Stokes equation, with a kinematic viscosity 



-1 



The second line in (28c) is a generalization of the convection term to systems not obeying 



Galilean invariance, where all combinations of V and factors second order in g transforming 
as vectors are allowed [H]. Finally, the last two lines describe couplings of the current den- 
sity g and gradients thereof to density gradients. Note that the density gradients in these 



coupling terms are all of the same generic structure (29). 



As already noted, the canonical model equations (28a) - (28c) conserve the total number 
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of particles. To make this explicit, we define 



p = pc + Ps 

V = Pc-ps 



(31a) 
(31b) 



where p denotes the overall particle density, and t] measures density difference between the 
two particle classes. The canonical model equations then attain the following form: 



dtp 
dtV 
dtg 



-Vg, 

-V ■g + p^ - {p + l)r]- p, 
-z/ig g g- -V{p + r]) 

1^2 4 



+ 



+ 



p_ 



•V)f 



d\pA) 



(V 

1 



1, 



V^g 



(32a) 
(32b) 
(32c) 



''d\pA 



[(V-g)a[p,r/]- (Vg + Vg*) d[pA]- 



The equation governing p expresses the overal conservation of particle number, whereas the 



source terms of equations (28a) and (28b) combine to determine the local dynamics of the 



relative density r] in (32b) 



Now we turn to the grand canonical model, where the single particle phase is coupled 
to a particle reservoir, resulting in a situation where single particles constitute an isotropic 
sea of particles which is maintained at a constant density p^. Particle number conservation 
is now violated, and the only non-trivial density dynamics takes place within the phase of 
cluster particles. The hydrodynamic equations corresponding to the grand canonical model 



can be obtained immediately by setting in (28a) - (28c) the density of single particles to a 
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constant value, yielding: 



ps 
dtPc 
dtg 



const. 



-V ■ g - Pe + (P° + Pc) P°, 



1^2 



-ypc 



V^g 



(33a) 
(33b) 
(33c) 



+ -(g-V)g+^ 

V2 ^2 



(V- 



pdpV2 



■Vpc 



Vp, 



4ul 



[(V ■ g)Vp, - (Vg + Vg*) Vpe] . 



One final remark is in order: The rescaling scheme introduced in equations (27a) - (27e) 



renders both, the canonical and grand canonical model equations virtually independent of 
particle shape. While these equations exhibit a weak dependence on the particles' aspect 
ratio L/d (via the rescaled collision integrals In,k), this dependence introduces only minor 
quantitative effects, which are negligible for all present purposes. To a good approximation 
we can thus set L/d = 1 while working with dimensionless variables, and assess the effects 
entailed by particle shape by restoring original units. Within our present approach, the 
effects of particle shape are purely quantitative, causing a numerical shift in the characteristic 
scales, but leaving the qualitative features of the problem unaffected. Deep within the 
ordered phase, i.e. for large densities, we indeed find a qualitative change of the ensuing 
hydrodynamic instability, as detailed in|V| Nevertheless, this statement has to be taken with 
a grain of salt because corresponding threshold densities are far beyond the validity of the 
hydrodynamic equations. 



IV. SPATIALLY HOMOGENEOUS SYSTEMS 



To investigate the implications of the hydrodynamic equations, we start with the simplest 
case by analyzing spatially homogeneous solutions. These considerations will provide the 
basis for the study of spatially inhomogeneous systems, which will be the subject of section 
|Vj Dropping all gradients, the hydrodynamic equations for spatially homogeneous systems 
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for the canonical model read 



dtp = 
dtV = 
dtg = 

For the grand canonical model we get 

dtpc = 
dtg = 



0, 

- (P + 1)^ - P, 
-z/ig g g. 

1^2 



(34a) 
(34b) 
(34c) 



-Pc + (p° + Pc) pI 

^2 



(35a) 
(35b) 



In both cases, the density dynamics decouples from the momentum current dynamics and 
can be addressed separately. 

In this section, our focus is on the stationary properties of the canonical and grand 
canonical model, respectively. While the dynamical approach to the stationary state is 
model dependent, the system's composition in terms of single particles and cluster particles, 



for given total density p, in the limit t — t- oo is identical in both cases [refer to (28a) - 



(28b) and (33a) - (33b)]. Since, moreover, the momentum current densities g obey identical 



dynamical equations, the ensuing analysis of the stationary state is equal for both models. 



A. Crossover to clustering 



To assess the density difference between the cluster particle and the single particle phase 



?7, we calculate the dynamical fixed point rj* of (34b), attracting the dynamics of ri(t) in the 
long time limit t — > oo: 



v*ip) 



P -P 
P + 1' 



(36) 



The defining equations (31a) - (31b) can be used to determine the corresponding (station- 
ary) fixed point densities of single particles [p*) and cluster particles [p*) as a function of 
the total density p. figure [3] summarizes these findings: Upon increasing the total density 
p, the ratio ri*/p continuously grows from ri*/p = —1 at p = 0, asymptotically approaching 
ri*/p = 1 as p — )■ oo. Based on the sign of t]*, two density regimes can be distinguished: 
In the low density regime (p <^ 1, r^* < 0) particle collisions, underlying the formation of 
clusters, occur at much smaller rates than cluster evaporation events. Only a small fraction 
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FIG. 3: Fixed points of the homogeneous equations for the grand canonical and canonical model: 
The stationary relative density rj* / as well as the stationary cluster and single particle density, 
and p*s/p, respectively. The larger the more cluster particles exist in the system. The 
vertical line corresponds to the density p above which the number of cluster particles exceeds the 
number of single particles. Note that < 1 holds for all finite values of the total particle density 
p. (This is of particular relevance in the context of the grand canonical model, where ps is can be 
considered as control parameter.) 



of all particles organize themselves in clusters leading to a relatively dense population of 
single particles and correspondingly small density of cluster particles. In the high density 
regime (p ^ 1, ?7* > 0), the situation is reversed: Large overall densities imply frequent 
particle collisions and, consequently, cluster formation and cluster growth dominate over 
cluster evaporation. In this regime, the number of cluster particles exceeds the number of 
single particles. 

The crossover between the single particle dominated low density regime and the cluster 
particle dominated high density regime occurs at the crossover density p = pb = 1, where 
both, the single particle and the cluster particle populations are of equal size [i.e. ?7*(p) = 
0]. The relation between the crossover density to clustering, and the geometrical shape of 
the constituent particles has been addressed previously in Ref. ^7\, based on agent-based 
simulations and a mean-field type analytical analysis. Using our definition of the crossover 
density p we can establish the corresponding relation simply by restoring original units 



[equation (26)]. Using packing fraction p ^ pLd instead of particle density, and assuming 
for the sake of simplicity L/d ^ 1, which allows us to estimate the particle surface Aq — Ld, 
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we find: 

eLd Sn eL . 

which correctly reproduces the findings of Ref. [37] (taking into account that the cluster 
evaporation rate is assumed to be proportional to the inverse particle length, e oc L~^). For 
the sake of completeness, we note that the definition of the clustering crossover density in 
reference is based on the cluster size distribution, and thus does not necessarily coincide 
with our definition. We stress, however, that in our description the scaling structure in 



equation (37) is completely generic. It is an immediate consequence of the characteristic 



scales of our model and of the fact that the rescaled hydrodynamic model equations are 



(virtually) independent of particle shape. The structure of equation (37) is thus robust 
under an arbitrary redefinition of the (rescaled) crossover density p. 

B. Homogeneous equations for momentum current density 

Having examined the composition of the system in terms of single particle and cluster 
particle densities, we now turn to a discussion of the spatially homogeneous solutions for 



the momentum current density g. Due to rotational invariance of (34c), only the magnitude 
g = |g| of the momentum current density, but not its direction, evolves in time. We can 
thus concentrate on the scalar equation 

d,g = -u,g-^g^ (38) 

which leads to the following fixed points g* as the attractor of the dynamics of g in the limit 
of long times: 

{0 forui > 0, 

, (39) 
90 = .^ foru,<0. 



It can be shown, that the coefficient in front of the cubic term in (38) is indeed strictly 



positive for all control parameters of density p and noise a consistent with ui < 0, ensuring 



the existence of the non-trivial fixed point in the second line of (39). 

Depending on the sign of the linear coefficient ui, two parameter regimes can thus be 
distinguished: Parameters leading to z^i > render stable an overall homogeneous and 
isotropic state with vanishing macroscopic fiow g = 0. Upon crossing the phase boundary 

z/i(p,a) = (40) 
20 



in parameter space, the isotropic solution gets unstable and a macroscopic current density 



of non-zero amplitude builds up. In equation (40) we used that the density difference rj, in 



the stationary limit, is a function of the total density p; cf. equation (36). Hence, in the 
limit of long times, ui is a function of the total density p and the noise parameter a, only. 



Using the definition of the coefficient z/i, equation (30a), we can readily calculate the 
shape of the phase boundary in the ap -plane: 



ae(p) = y-2 In + p-2j , (p > v^), (41) 
where we used /i_o = — | and Ji^i = ^. The corresponding phase diagram is shown in figure 

a 

To conclude this section, we note that the analysis of spatially homogeneous systems 
corroborates the general physical picture of active systems, that was alluded to in the in- 
troduction (e.g. cf. reference PT]): Even in the absence of noise, a = 0, for which the 
threshold density p^'^^ is lowest, the fully isotropic state g = remains stable up to a critical 
density p^^^ {a = 0) = y/S p, which lies well beyond the density p indicating the crossover to 
clustering. We thus extract the following physical picture; cf. figure |4j For low densities, 
p < p, cluster evaporation dominates over cluster assembly via particle collisions and clusters 
form only transiently. The system most closely resembles a structureless, isotropic "sea of 
particles". At intermediate densities, p < p < p'^^\ particle collisions are more frequent. The 
emergence of clusters is now a virtually persistent phenomenon, with cluster evaporation oc- 
curring at a lower rate than cluster formation and growth. Yet, the collision rates between 
clusters (i.e. collisions among cluster particles) are still too low to orchestrate macroscopic 
order, leading to an overall isotropic "sea of clusters". Finally, for large densities, p > p'^^\ 
the frequency of collisions among clusters is high enough to establish collective motion even 
on macroscopic scales. 



V. STABILITY OF INHOMOGENEOUS HYDRODYNAMIC EQUATIONS 

From our hitherto discussions, we have ascertained that the isotropic, homogeneous state 
(p = const, and g = 0) becomes unstable for sufficiently large densities. Yet, from a 
purely homogeneous analysis we cannot tell anything about the spatial structure of such a 
macroscopic broken-symmetry state. Nor can we be sure that the isotropic and homogeneous 
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FIG. 4: Phase diagram given by the homogeneous equations for the canonical and grand canonical 
model. For noise values smaller than the critical value, a < cTc, the isotropic state becomes unstable, 
giving rise to a state of collective motion of non-zero macroscopic momentum current. For a > Uc 
the isotropic state {qq = 0) represents a stable solution. The vertical dotted line indicates the 
transition density p'^^^ at zero collision noise u = 0, and the vertical dashed line corresponds to 
the crossover density p, above which the number of cluster particles exceeds the number of single 
particles. 

solution for p < p'^{cr) is indeed stable with respect to spatially inhomogeneous perturbations. 
In this section, we therefore test the linear stability of the homogeneous isotropic and non- 
isotropic base states with respect to wavelike perturbations of arbitrary wave number. Unlike 
the homogeneous model equations, the full hydrodynamic model equations are different 
for both, the canonical and grand canonical model, implying different dispersion relations 
describing the growth of such wave-like perturbations. We will thus analyze both models 
separately, and show that particle conservation does indeed influence pattern formation in 
essential respects. 

A. Linearization about stationary, spatially homogeneous base states 

We start by linearizing the hydrodynamic equations for the canonical model. In the 
canonical model, the total number of particles is conserved, and the appropriate base state 
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reads (cf. section IV) 



P 
V 



Ph. = const. 



Ph- Ph 



Ph + l 



V1U2 



"9' 



(42a) 
(42b) 

(42c) 



where denotes the unit vector in the direction of the homogeneous polarization, and 
where all fields of the base states are assumed to be constant both in space and time. We 



are going to investigate the linear stability of the solutions (42a) - (42c) against wave-like 
perturbations, employing the following ansatz: 



p(x,t) = p/, + (5p(x,t), 
r/(x,t) = r]* + 5r/(x,t), 
g(x,t) = gft, + (5g(x,t). 



(43a) 
(43b) 
(43c) 



where the perturbations are plane waves 



(44a) 
(44b) 
(44c) 

In the equations above, q denotes the wave vector and s is the growth rate. Inserting this 



5p(x, 


t) 


= Spoe^'- 


h«q-x 


(5r7(x, 


t) 


= Sme''' 


-jqx 


5g(x, 


t) 


= 5goe''- 


f jq-x 



ansatz into the hydrodynamic equations (32a) - (32c), we obtain the following eigenvalue 
problem: 



sSpo 
s5r]o 
s5go 



{2ph -1]* -1) 5po - (1 + ph) 5r]o - iq ■ 5go 

2 

2 



dpT^2P (2 , / ■ N 2 \ ct 

' "■ "h Sh + [gh -i^Jgh- ^gh] - OpVi gh - 



4 



5p 



(45a) 
(45b) 
(45c) 



i^ShSh + [gh ■ ^q) g/. - ^ gh ) - dr^^i Sh - 



zq 

T 



61] 



/■ 2 

(,+ N q 2 

— «q -Sh) - gh 

1/2 4z/2 Z/2 



X 2pfi: 

Sgo gh [gh ■ dgo) 

^2 



+ 



1^2 



gh{i(\ ■ 5go) - i(\{gh ■ 5ga 
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Unlike the canonical model, the grand canonical model conserves the number of single 
particles, but not the total number of particles. The appropriate base state in this case 
reads: 



Ps 
Pc 



const 

P*ciPs. 



Pt 



Z/ii/2 



(46a) 
(46b) 

(46c) 



We investigate the linear stability of these solutions, using a perturbation ansatz analogous 



to equations. (43a) - (44c) 



Pc(x,t) = p* + (5pc(x,t). 



with 



st+iq-x 



(5pc(x,t) = (5p^e 



Inserting this ansatz into equations (33a) - (33c), we obtain: 



sSpl 



(p, - 1) 5p° - iq ■ Sgo 

dp,U2P / 2 , / ■ N 2 
' f^ShSh + (gh ■ ^q) gh - ^gh 



+ 



c 



jJiKj 2 

— i«q ■ Sh) - - gh 

Z/2 4z/2 Z/2 



(47a) 
(47b) 



(48a) 
(48b) 



(49a) 
(49b) 



5go - — g/i (g/i • 5go) 



1^2 



Z/2 L 



^^(iq ■ Sgo) - iq(g/i • 6go) 



B. Stability of the disordered state go = 

We start by considering the homogeneous and isotropic base state, which was shown to 
be stable against spatially homogeneous perturbations for p < p^'^\a); cf. section 



IV B 



To assess the stability of this state with respect to perturbations of arbitrary (non-zero) 
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FIG. 5: Fastest growth rate of 5R[s(g)] as a function of of the wave number q = |q| for the canonical 
(a) and grand canonical model (b), each for o" = 0. The disordered state is stable for all wave 
numbers q if p < p^'^\cr). The marginal case, p = p^'^\a), is dashed. An instability (grey) occurs 
for densities larger than the corresponding homogeneous critical density p^^\a). Similar behavior 
is found for £T / 0. 



wavevectors in the canonical model, we use the linearized hydrodynamic equations (45a) 



- (|45c|) with Qh = 0. The resulting eigenvalue problem is most conveniently expressed in 



matrix form: 



\S9o J 





2p,,-r/*-l -{1 + Ph) 
—iq/A —iq/A 



—iq 
—iq 

J2 



z/i - g7(4z/2) j 



(50) 



The corresponding eigenvalue problem for the grand canonical model is found from equations 

-S i. 



(49a) - (49b), and attains the following form: 



r 



—iq 



-zg/2 -z^i - gV(4z/2) 




(51) 



For 



0, (45c) or (49b), respectively, implies q||5go, allowing to replace the vectors q 



and 5go by their respective magnitudes q and dg^. We solved both eigenvalue problems 
numerically, for arbitrary wavenumbers g > with the results shown in figure [5j Note that 
in both models the real parts of all eigenvalues are negative for all wavenumbers g > 0, 
provided the particle density p < p^^\ The spatially homogeneous, isotropic state is thus 
stable against small perturbations with arbitrary wavevectors. 

For densities p > p^^\ in contrast, a narrow band of positive eigenvalues emerges in both 



models, located at wavenumbers g <^ 1. Equations (50) and (51) evaluated at g = return 
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nothing but the hnearized versions of the homogeneous hydrodynamic equations, (34a) 



(34c) and (35a) - (35b). To gain new insights, we will therefore examine the limit g — )■ and 
consider the eigenvalues of the above coefficient matrices to leading order in the wavenumber 
9- 



In this limit of small wavenumbers, the grand canonical coefficient matrix, given in (51), 
approaches diagonal form and the dynamics of density fluctuations 6p^ and momentum cur- 



rent density fluctuations 6 go practically decouple. Since < 1 (cf. flgure [3]), the flrst 
eigenvalue s^^'^^ = ps — 1 + C(g^) is strictly negative and density fluctuations decay exponen- 
tially. The second eigenvalue, s'^'^^ = —ui + 0{q'^), is positive at small wavenumbers leading 
to an instability in the momentum current density against long wavelength fluctuations. 



In the case of the canonical model (50), the coefficient matrix approaches block diag- 
onal form in the limit of small wavenumbers. Again, the dynamics of momentum cur- 
rent density fluctuations 6go practically decouples from density fluctuations (5po and Srjo), 
with momentum current density fluctuations being amplifled by virtue a positive eigenvalue 
Sg*^'* = — z/i + at small wavenumbers. In contrast to the grand canonical model, how- 

ever, particle conservation entails a marginally stable mode s^^\q = 0) = 0, which turns 
positive for g > 0: s^^^ oc (the remaining eigenvalue ^2*^^ = —(1 + ph) + C^(g^) is strictly 
negative). 

To sum up, the study of the linear stability of the homogeneous, isotropic state against 
spatially inhomogeneous perturbations of arbitrary wave vectors strongly suggests that par- 
ticle conservation plays a vital role in the context of pattern formation. Both models exhibit 
spontaneous symmetry breaking by establishing a state of macroscopic collective motion. In 
the canonical model, in addition, conservation of total particle number entails a marginally 
stable density mode at g = which is absent in the grand canonical model. This mode, in 
turn, gives rise to a density instability at small, non-zero wavenumbers, accompanying the 
spontaneous symmetry breaking event for p > p^^\ We note, however, that, at this point 
of the discussions, the existence of a narrow band of unstable modes at small wavenumbers 
does not allow for any conclusions concerning the structure of the macroscopic density and 
momentum current density for p > p^^\ We will address this issue in greater detail in the 
following section. 
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C. Stability of the broken symmetry state go > 



Both, the canonical and grand canonical model exhibit spontaneous symmetry breaking 
for overall densities p > p^^^a). To illuminate the spatial structure of this broken symme- 
try state, we start from the most simple case of a spatially homogeneous state of collective 
motion, and examine its stability with respect to wavelike perturbations in the hydrody- 
namic particle and momentum current densities. Without loss of generality, we assume the 
direction of the macroscopic momentum current density to coincide with the x-direction 
and choose = go e^.. The wave vector q of the perturbation fields is assumed to make 
an angle ip with the macroscopic momentum current density gh, yielding ip = Z(q, e^) and 
q = g (cos ('?/'), sin (ip)) with q = |q|. 



The linearized canonical model equations (45a) - (45c) then attain the following form 



s6po = -iqcos{ilj)6gj:fi -iqsm{ij;)Sgyfi, 
s6r]o = {2ph - 11* (ph) - 1) 5po - (1 + Po) % 
-iqcos{i!)5g.^^o - iq^^^W^gyfi, 

sSg^fi = 



iqcosiip) ( gl ^dpU2 - ^ ) - ( dpUi 



+ 



-iq cos{i)) [gl ^d^U2 - 

+ ^^^^°^(^)^^°-4^ 
+zq—sm{ilj)goSgyfi, 

^2 



2pK 
1^2 



x,0 



s5gyfl = -^igsin(V^) (^gl ^dp^2 + ^ ^Po 

-^igsin(V^) (^l ^dv^2 + ]^ 5t]q 

(_ 

-iqsm{tp)— go 6g^^o 

^2 

+ ( zgcos(V^)— - ) bgy^o- 
V V2 4i/2 / 



—9oOpV2 



- 9o ( dr^^i - ^9ld7j^2 



5po 
5r]o 



(52a) 
(52b) 

(52c) 



(52d) 
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The corresponding equations for the grand canonical model read: 



{ps - l)5pl - iqcos{ij)5gxfi - iqsm{^)5gyfi, 
Uqcosi'^) (gl ^dph 



1 



(53a) 
(53b) 



- ^0 ( dp^vi - ^gldp^U2 



5pl 



+ [ iq cosOj)— go 

^2 



4:U2 



-9o (^9x,Q 



V2 



+iq— sin(^)5(o dgy^, 

V2 



sSgyfl = -^igsin(V') (^go ^dpV2 + 1^ 5p 

(_ 

-iq sin^ip)— go 6g^^o 

^2 

+ Lcos(^)^^7o-|;^J 5gy,o- 



(53c) 



In equations (52a) - (53c) we used —vi 



^ ^2 
I/O yo 



0, which directly follows from the defini- 



tion of go, given in equation (39). We numerically solved both eigenvalue problems in the 
immediate vicinity of the ordering transition line 

In the case of the canonical model, we find that the most unstable mode occurs for 
longitudinal perturbations, i.e. perturbations with wave vectors parallel to the direction of 
macroscopic motion, q || go (V' = 0). |6^ shows the corresponding eigenvalues as functions of 
the wavenumber q for a set of density values slightly beyond p = p^^\ Further inspection of 



the coupling coefficients in equations (52a) - (52d) reveals that this longitudinal instability 



only affects the amplitude of g leaving the direction unchanged: For ip = 0, the dynamics of 
Sgyfi decouples and momentum current density fluctuations perpendicular to the direction 
of macroscopic motion decay exponentially. 



with a rate 



3f? 



.(c) 



< 0, 



(54) 



(55) 



which approaches zero for g — )■ 0, as expected for a broken symmetry variable. To assess 
the nature of the instability in greater detail, we calculated the eigenvector corresponding 
to the most unstable longitudinal mode (evaluated at the most unstable wavenumber). 
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FIG. 6: Largest growth rate of 3^[,s(g)] as a function of the wave number q for o" = and several 
values for the total particle density p. The marginal p = p^^^ is dashed. Further values are p = p^'^^ + 
A with A indicated in the figure, (a) Canonical model: For p > p^^\ longitudinal perturbations 
{ip = 0) are unstable, (b) Grand canonical model: Transversal {{ipl = it/ 2) perturbations are 
unstable closely above the critical density p*^*^-* (refer to green and blue curve corresponding to 
A G {0.05,0.5}). For larger densities, i.e. A > 0.7 for cr = 0, the transversal instability re- 
stabilizes again (dotted curves, A G {0.8, 1}). However, the density regime hosting this transversal 
instability vanishes completely for noise values larger than ar, as illustrated in[7j 

It turns out that this eigenvector has approximately equally large components along the 
remaining fluctuation amplitudes Sg^fl, Spo and 6riQ. This is consistent with our previous 



findings, indicating that the density mode, which was alluded to in section VB and which 
turns unstable at p = pc, renders the state of homogeneous collective motion unstable to 
fluctuations of the magnitude of the momentum current density. We further note that this 
picture is in agreement with previous numerical [3T] and analytical [22] results (cf. [T]). 

The stability regions of the grand canonical model strongly deviate from the above pic- 
ture. Setting ip = (longitudinal perturbations), we calculated the largest eigenvalue s^°^ 



of the linear system of equations (53a) - (53c) 



4[(l-p,)2 + p2(i| + |e--^)]' ^^^J 

which is always negative since p^ < 1. In contrast to the canonical model, longitudinal per- 
turbations thus always decay exponentially fast in the grand canonical model. For perturba- 
tions in transverse directions, in contrast, a positive eigenvalue can be found for sufficiently 
low noise levels a < a^, with the fastest growing modes posessing wavevectors q ± go- [S)? 
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shows the eigenvalue of the most unstable modes, which occur for o" = 0. To assess the 
implications of this instability for the dynamics of the various fluctuation amplitudes, we 
numerically examined the eigenvector corresponding to the positive eigenvalue, evaluated at 
the most unstable wavenumber. For densities in the vicinity of the ordering transition, we 
find that this eigenvector has approximately equal components in both momentum current 
density fluctuation amplitudes, Sgxfl and Sgy^o, but an essentially vanishing component along 
the direction of density fluctuations 6p^. The corresponding instability can thus be classified 
as a hybrid shear/splay instability, leaving the spatially homogeneously distributed particle 
density virtually unaffected. 

Three remarks are in order: First of all, for noise values ac{p — oo) > a > ar, the state 
of homogeneous collective motion becomes linearly stable with respect to arbitrary pertur- 
bations, including transverse perturbations. Secondly, even the most unstable eigenvalues 
"restabilize" for densities which are in the vicinity of the ordering transition threshold, which 
is depicted in [7j Finally, a restabilization can also be "observed" for the longitudinal insta- 
bility in the canonical model. In this case, however, the restabilization occurs for relatively 



large densities and thus lies outside the range of validity of the linearized equations (52a) 
(pdl). 



VI. DISCUSSION AND CONCLUSION 



To conclude, we discuss and summarize our main findings. To study the onset of collective 
motion in active media, we started out with a simplified model for a system of self-propelled 
rod-like particles of variable aspect ratio. Collective motion was assumed to be established 
in a completely self-organized fashion solely by means of interactions among the constituent 
particles and in the absence of any external alignment fields. These interactions were as- 
sumed to occur via binary, inelastic particle collisions during which the rods align their 
direction of motion. Moreover, interactions were assumed to be subject to noise, which we 
controlled by a single model parameter a. To assess some of the structural properties of 
such systems, we associated each of the particles with one of two classes: single particles and 
cluster particles, each with the corresponding density fields denoted as ps and pc- The class 
of cluster particles hosts all particles belonging to some coherently moving group of particles, 
which we referred to as cluster. The rest of the particles can be imagined to make up an 
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FIG. 7: Phase diagram determined from the homogeneous equations of the grand canonical model, 
as a function of noise level a and single particle density p^, now complemented by the results 
obtained from the stability analysis of the linearized inhomogeneous equations: Whereas longitu- 
dinal perturbations decay within the homogeneous phase boundary, there is a zone (grey shaded) 
where transversal modes become linearly unstable. The width of this zone gradually decreases for 
increasing noise values a, and vanishes completely above some critical noise value dr (horizontally 
dotted line). 

isotropic sea of particles and are associated with the class of single particles. Using this clas- 
sification scheme, we implemented simple interaction rules, representing cluster nucleation, 
cluster growth and cluster evaporation; the latter is assumed to occur at some fixed rate e. 

To illuminate the self-organization of collective motion, we set up an analytical, kinetic 
description of such systems, focusing on two archetypical modeling frameworks. Firstly, we 
considered isolated systems in which the total number of constituent particles is a conserved 
quantity. This case was referred to as the canonical model. Secondly, we examined open 
systems, which we referred to as the grand canonical model. Open systems are in contact 
with a particle reservoir which keeps the density of single particles at a constant level. 

Inspecting the corresponding hydrodynamic equations, we were able to establish the fol- 
lowing physical picture, portraying the formation of collective motion via dissipative particle 
interactions: For both, the canonical and the grand canonical model, we identified two char- 
acteristic density scales p and p^'^^{a), with p^'^\ci) > p, which allowed us to distinguish three 
density regimes. 

For low densities, p < p, the rate at which particles collide is much smaller than the rate 
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at which clusters disassemble. In terms of a particle based picture, this regime corresponds 
to a situation, where particle clusters are unstable, evaporating shortly after their nucleation. 
In the stationary state, the vast majority of particles populates the single particle phase, 
rendering the system homogeneous and isotropic even on mesosopic scales. This low density 
regime terminates at the characteristic density p, where both classes exchange particles at 
equal rates. 

In the contiguous regime of intermediate densities, p < p < p^^\ the overall rate of cluster 
formation and growth outstrips the rate at which clusters evaporate, and the majority of 
particles becomes organized in clusters. Translated to a particle based notion, clusters 
grow to finite sizes and persist over macroscopic time scales. Clusters of coherently moving 
particles now dominate the physical picture on mesoscopic scales. Yet, interactions among 
clusters are too rare to establish a macroscopic state of collective motion. On hydrodynamic 
length scales, the system can be viewed as a homogeneous and isotropic sea of clusters. 

For densities exceeding the critical density, p > p^^-'(cr), collisions within the cluster 
phase occur at sufficiently high rates, and macroscopic collective motion emerges. The 
homogeneous and isotropic state, which has been shown to be stable within the two preceding 
regimes, thus gets unstable and rotational symmetry is spontaneously broken. Resorting 
to a particle based image, we can imagine the mean cluster size to reach a "percolation 
threshold" , leading to coagulation and net alignment between clusters. 

While the qualitative features of the canonical and the grand canonical model are the 
same in the low and the intermediate density regime, the establishment of collective mo- 
tion in the high density regime differs in important respects in both models. We found 
that in the grand canonical model, a broadly extended region in parameter space exists, 
where a spatially homogeneous state of macroscopic collective motion exists and is actually 
stable. Except density, the key parameter controlling the stability of a spatially homoge- 
neous flowing state is the noise amplitude a. For low noise levels the homogeneous flowing 
state gets unstable toward transverse perturbations (i.e. perturbations with wavevectors 
q perpendicular to the direction of the macroscopic flow). We note, however, that these 
instabilities are remarkably weak, i.e. the corresponding growth rates are smaller than those 
of the longitudinal instability by a factor of ~ 10 (cf. flgure[7]), and "restabilization" of 
the spatially homogeneous flowing state occurs upon increasing the density only slightly 
beyond the threshold p^^\a). Interestingly, this transverse instability vanishes altogether, if 
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angular diffusion is slightly enhanced upon increasing a. Hence, for intermediate values of 
cr, the system directly establishes a homogeneous state of collective motion, which is stable 
against arbitrary perturbations of small magnitude. Finally, if the noise is too strong, order 
is destroyed and the system remains isotropic even for arbitrarily large densities. This last 
statement is, of course, shared among all active systems [29], particle conserving or not, and 
thus applies equally well to the canonical model. 

In the case of the canonical model, a spatially homogeneous base state is unstable toward 
longitudinal perturbations (i.e. perturbations with wavevectors q parallel to the direction 
of the macroscopic flow) for all values of the noise parameter a. Both, the magnitude of 
the macroscopic velocity field and the particle density are prone to this kind of instability. 
This is in agreement with previous analytical [22] and numerical [31] results for particle 
conserving systems, where the emergence of solitary wave structures has been reported in 
the vicinity of the ordering transition p > p^^\a). The longitudinal instability thus seems 
to be a quite generic feature of particle conserving systems with hard core interactions. For 
an interesting counter example we refer the reader to Ref . [SB] , where a particle conserving 
system with topological interactions has been studied. 

We can now combine our findings for both, the canonical and the grand canonical model, 
to offer the following mechanistic explanation concerning the emergence of the longitudinal 
instability. The prerequisite, underlying the establishment of coherent motion, is embodied 
by two basic processes: Cluster nucleation by collisions among single particles, and cluster 
growth by alignment of single particles to clusters. Only if, by virtue of these processes, the 
concentration of cluster particles grows sufficiently large, clusters are able to synchronize 
their movements by coagulation and macroscopic collective motion emerges. 

Now consider the effect of a density fluctuation in an otherwise homogeneous state of 
macroscopic collective motion. In the grand canonical model, where the density of single 
particles is kept fixed by virtue of a particle reservoir, this fluctuation occurs within the 



class of cluster particles. We can use the right hand side of equation (35a), to assess the 
implications of such a fluctuation on the local composition of the system in terms of cluster 
particles and single particles: 

(Ps + Pc) P° = Pc (57) 
Note that this equation captures the balance of the two particle currents between the single 
particle and the cluster particle phase in the stationary limit. As can be seen from this equa- 
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tion, locally enhancing the density of cluster particles implies a net current from the cluster 
particle phase into the single particle phase, thus counteracting the effect of the original 
density fluctuation. Conversely, locally diminishing the density of cluster particles leads to 
the opposite effect. Density fluctuations are thus damped in the grand canonical model and 
do not impact the macroscopic velocity field, which is set up by the cluster particles. 

Exactly the opposite happens in the particle conserving canonical model. Again, consider 
a spatially homogeneous base state of macroscopic collective motion. Particles are then 
distributed among the phases of cluster particles and single particles as determined by the 
balance equation [cf. ( |34b ) 

p[p-n) = p + r], (58) 

where the left hand side describes cluster nucleation and condensation, and the right hand 
side corresponds to cluster evaporation. This can be seen by using the definitions of the 
relative density rj = pc — Ps, and the total particle density p = ps + Pc- Now, consider a 
fluctuation in the total density p, where, for the sake of simplicity, we assume the relative 
density r] to remain constant. In regions, where the fluctuation leads to an increase in the 
total density by a factor > 1 we have 

kp{kp — rj) > kp + rj. (59) 

Hence, the particle current into the cluster particle phase grows. As a consequence, the 
local value of the momentum current density increases, since the cluster particles are the 
"carriers" of the macroscopic momentum. In contrast, in regions, where the fluctuation 
decreases the total density by a factor k' < 1 we have 

k'p{k'p -T]) < k'p + T]. (60) 

There the cluster particle phase gets depleted and the local magnitude of the momentum 
current density declines. As a result, high density regions move at faster speeds than low 
density regions, gathering more and more particles on their way through the system. Con- 
versely, lower density regions continually lose particles to the faster high density structures. 
In particle conserving systems, every density fluctuation thus automatically triggers a corre- 
sponding fluctuation in the momentum current density, which in turn amplifies the density 
fluctuation. As a result of this process, high density bands of collectively moving cluster 
particles might emerge [31j. These bands being interspersed by regions where the particle 
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density has fallen below the critical density p^'^^ (and possibly below p), leading to local 
destruction of clusters and collective motion. 

We close by adding some remarks on the importance of the particles' shape on the estab- 
lishment of collective motion on hydrodynamic scales. We found that the impact of particle 
shape on the macroscopic properties of such systems is purely quantitative in the framework 
of our present study: Varying the particles' aspect ratio results in a shift of the character- 



istic density scales p and p^^\a), which we quantified in equation (37). Qualitatively, our 
conclusions concerning the macroscopic properties of these systems remain unaffected by a 
change in the particles' aspect ratio. Note that, in our approach, the aspect ratio basically 
determines the total scattering cross section and thus "merely" impacts the rate at which 
particles collide. We stress, however, that in real systems particle shape is likely to have a 
profound impact on the entire physical picture of particle interactions, and not just on their 
rate. The study of those effects lies outside the scope of our present work and would be an 
interesting topic for future research. 
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Appendix A: Derivation of Boltzmann collision cylinder for driven rods 

In the framework of our Boltzmann-like description, binary collisions, such as equations 
(|3|, ^ and ([T]), occur with a certain rate T depending on particle shape (L and d), relative 
angle of both collision partners 612 = \0i — 02\ and the constant velocity v. The quantity 
r(L, d, 612) characterizes the collision area per unit time - more commonly referred to as 
Boltzmann collision cylinder. On the scale of the Boltzmann equation, binary collisions 
occur locally, say in an infinitesimal volume element centered at r. Assume that particle 1 
has an orientation 61. Then, r{L,d,6i2) dt gives the area around particle 1 in which every 
particle with orientation 62 will collide during a time interval [t, t + dt] with particle 1 . As 
a consequence T{L, d, Ou) /(r, Oi, t) /(r, 6*2, t)d0id62 equals the number of collisions per unit 
time and unit area at time t, with f{r,6,t) denoting the one-particle distribution function. 

To determine r(L, d, 612), we take a microscopic point of view. Since the model employed 
in this work assigns to each particle a velocity vector pointing along its rod axis, we can 
distinguish "head" and "tail" . Referring to figure [8} without loss of generality we assume 
n — 9i2 = E [0, tt] (negative relative angles lead to the same result), and consider the blue 




FIG. 8: Illustration of the collision cylinder in the rest frame of the blue rod. The red lines 
indicate the excluded volume due to the finite expansion of the rods. THe quantity v^-ei denotes 
the magnitude of the relative velocity of those rods making a relative angle Qyi = it — 9 with the 
blue rod's axis, and is given by Vrei = v\e{9) — e(0)| = 2v\ sin(0i2/2)|. 
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rod, with the position of its head indicated by the blue dot. All rods of relative orientation 
6'i2 = ^1 — ^2) and with their heads lying in the area S* = A U U 5*2 at time t, will collide 
with the blue rod during the time interval [t,t + dt]. Since A, Si and 5*2 are disjoint, 



\S\ = \B\ + \Si\ 



where |X| denotes the area of the region X. The respective areas are given by: 



\A\ = dt Vrei {L — d) \ sin 6\ = dt Vrei {L — d) \ sin 6 



12 



and 



|5'2| + ISil = dtVrel d 



TT-e 



simd) + 



2 dt Vrpi d . 



(Al) 



(A2) 



(A3) 



Returning to the laboratory frame we have Vrei = ^|e(^^i) 
that T=\S\/dt (cf. eq. (lAl])), we find: 



6^2 



2w|sin(ei2/2)|. Noting 



T{L,d,9i2) =^vd 



sm 



712 

2 



1 



L/d 



\sin6 



12 



(A4) 



In figure [9| T(L,d,9i2) is shown as a function of relative angle 6^12 for different particle 
lengths, whereby the particle width d is kept fixed. Increasing L/d shifts the most probable 
collision from 612 = vr for L/d = 1 (the case of a sphere; 612 = tt leads to the largest value 
of the relative velocity), towards 612 = vr/2 for L/rf — )■ 00 (limiting case of a needle; largest 
target area for 9i2 = 7r/2). 



Appendix B: Derivation of the gradient terms in the hydrodynamic equations 

To assist the reader in tracing back the emergence of the gradient terms in the hydrody- 
namic equations (28a) - ( 28c[ ) [and, likewise, in equations (32a) - (32c) and (33a) - (33c)], 
we briefly summarize the main steps in the derivation of these equations. All gradient 
terms in the hydrodynamic equations ultimately arise from the convection term in the first 



line of equation (24c) and the closure relation obtained by quasi-statically approximating 



(24d). Here we collect all such (complex) gradient terms and give a brief derivation of their 



vector-analytic counterparts. As in the main text, we identify C and M^, i.e. 



/ = /.. + ^/v e C ^ f 




(Bl) 



To distinguish (genuinely) complex from purely real quantities, we assume / G C and p G 
in the following. 
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FIG. 9: T(L,d,6i2) as a function of the relative angle 9i2 for different values of aspect ratio 
For the figure, we chose for particle width d = 1 and for particle velocity v = 1. Increasing the 
aspect ratio L/d, the most probable collision approaches 612 = it/2, whereas for L/d = 1 the most 
probable collision is the head-head collision with 612 = vr. 

{dx + idy)p 



Using (Bl) we immediately obtain 



{dx + idy)p = Vp. 



(B2) 



{dx -idy){dx + idy)f 



By straightforward expansion we find 



{dx - ldy){dx + ldy)f = {dl + dl)f ^ VH. 



(B3) 



{dx - idy)f 

Decomposing / into real and imaginary part and expanding, we find 

{dx - ldy){fl - fl + 2lfxfy) = dxfl - dxf; + 2dy{fxfy) 

+^ [dyfy-dyfl + 2dx{fxfy 

1 

2 



= 2 



^^{f^f,) - i5.,5.f2 



2f (V ■ f) + 2(f ■ V)f - Vf^ 



(B4) 
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where denotes the j-th Cartesian unit vector. 

[{dx + idy)f] [{dx - idy)p] 

Expanding and collecting real and imaginary parts, we find 

[{dx + idy)f][{d:c-idy)p] = d^f^d^p - dyfyd^p + d^^ydyP + dyf^dyp 

+ i {-d^cfxdyP + dyfydyP + d^fyd^cP + dyfAp) . 

Thus, 

m + idy)f] m - idy)p] = {idifj)dip + idjfi)dip)ej - (V • f ) Vp (B5) 

= [(Vf) + (Vf)*]Vp-(V-f)Vp. 

p{dx-idy)p 
We find 

f{d, - ldy)p = fld^P - fld^P + 2UydyP + I {-f'^dyP + f^yP + 2fjyd,p) . 

Hence, 

fid, - idy)p = 2fjjdjpe, - f^Vp - 2 f (f • Vp) - f^Vp. (B6) 
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